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The many-body Coulomb repulsive energy of strictly correlated electrons provides direct informa- 
tion of the exact Hohenberg-Kohn exchange-correlation functional in the strong interaction limit. 
Until now the treatment of strictly correlated electrons is based on the calculation of co-motion 
functions with the help of semi-analytic formulations. This procedure is system specific and has 
been limited to spherically symmetric atoms and strictly ID systems. We develop a nested opti- 
mization method which solves the Kantorovich dual problem directly, and thus facilitates a general 
treatment of strictly correlated electrons for systems including atoms and small molecules. 
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I. INTRODUCTION 

Kohn-Sham density functional theory (KSDFT) is the 
most widely used electronic structure theory for systems 
with many electrons, from gas phase molecules to con- 
densed matter systems in particular^'^. It is in prin- 
ciple an exact theory, and is able to yield the exact 
ground state energy and density using a fictitious sys- 
tem of non-interacting electrons. The key component of 
KSDFT is the exchange-correlation functional. Tremen- 
dous progress has been made in the past three decades for 
constructing approximate exchange-correlation function- 
al based on the known information from the uniform 
electron gas"^^^. However, such approximate exchange- 
correlation functionals are known to fail for strongly cor- 
related systems, such as the chromium dimcr^ or Mott- 
Hubbard insulators^'^. Several recent studies indicate 
that the construction of exchange-correlation functionals 
for general strongly correlated systems can be extremely 
difhcuh^-". 

Recently, the behavior of the exchange-correlation 
functional has been revealed in the limit of strictly cor- 
related electrons (SCE)^^"^^. The many-body Coulomb 
repulsive energy of SCE determines the exact exchange- 
correlation functional in the strong interaction limit, 
without artificially breaking any symmetry of the system 
or introducing any tunable parameters. The information 
provided by the SCE limit is complementary to that pro- 
vided by the uniform electron gas. For given electron 
density profile, the SCE limit is described by minimizing 
the many-body Coulomb interaction energy with respect 
to all antisymmetric wavefunctions in a 3A^ dimensional 
space, where N is the number of electrons in the sys- 
tem, under the additional constraint that the wavefunc- 
tion is consistent with the electron density^^. Mathe- 
matically, this daunting minimization task is an optimal 
transport problem^^ with Coulomb cost function^'^'^"'''^^. 
From physical intuition, this problem can be simplified 
by introducing N co-motion functions : — >■ M'^^^. 
These co-motion functions characterize the relative posi- 



tions of all the electrons with respect to one given elec- 
tron in the SCE limit. To the extent of our knowledge, in 
practice the co-motion functions can only be determined 
for one dimensional systems^^ and spherically symmet- 
ric atoms^^'^'^'^^, with the help of semi-analytic meth- 
ods. Little is known about the shape or even the ex- 
istence of the co-motion functions for general systems 
including small molecules. On the other hand, the opti- 
mal transport problem with Coulomb cost function can 
be equivalently solved by its dual formulation, called 
the Kantorovich dual problem^-^'^^'^^"^". The main ad- 
vantage of the Kantorovich dual problem is its poten- 
tial applicability to general systems, ranging from atoms 
and molecules to condensed matter systems. However, 
the Kantorovich dual problem is formulated as a maxi- 
mization problem with an infinite number of constraints, 
which is impossible to be implemented directly. These 
limitations severely restrict the applicability of the SCE 
limit to systems of practical interest. In this paper, 
we develop a novel method that solves the Kantorovich 
dual problem directly. We overcome the difficulty of in- 
finite number of constraints via a nested optimization 
approach. Our method provides a more general treat- 
ment of the exchange-correlation functional in the SCE 
limit for atoms and small molecules. 



The rest of the paper is organized as follows: in Sec- 
tion II, we briefly review the SCE limit, the optimal 
transport formulation and the Kantorovich dual problem, 
and present the nested optimization method for solving 
the Kantorovich dual problem directly. In Section III, we 
establish the applicability and accuracy of this method 
for the 3D beryllium atom and a model quantum wire 
system in ID, for which accurate results can be obtained 
semi- analytically using the co- motion formulation. Next, 
we demonstrate the applicability of our method to a 
model trimcr with various number of electrons in 3D, 
for which the SCE limit cannot be calculated by existing 
techniques. The conclusion and future work is given in 
Section IV. 
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II. THEORY 



According to the Hohenberg-Kohn theorem^, the 
ground state energy of a system can be obtained by min- 
imizing the following functional with respect to the elec- 
tron density p(r): 



E[p] = F[p] 



t(r)p(r) dr. 



(1) 



Here Wcxt(r) is the external potential, and F[p\ is the in- 
ternal energy functional, which is a universal functional 
of the electron density and consists of the kinetic energy 
and the Coulomb repulsive energy between the electrons. 
Formally F[p\ is defined by minimizing over all the an- 
tisymmetric wavefunctions ^' which are consistent with 
p{y) as 



F[p] 



mm 



(2) 



Here T — — Si=i kinetic energy operator and 

Vee = Y.iLiY.^>i\^i ~ i^iT^ is the Coulomb repulsive 
energy operator. A,; is the Laplacian operator on the i- 
th electron. The strong interaction limit considers the 
situation when the Coulomb repulsive energy dominates 
over the kinetic energy, in which case the internal energy 
functional can be approximated as^'^ 



F[p] 



mm 

TJp] 



* T 



mm 



(3) 



The first term Ts[p] is the Kohn-Shani (KS) kinetic en- 
ergy functional corresponding to a non-interacting in- 
dependent particle system^. The second term Vgg^^[p] 
is the minimal Coulomb repulsive energy among all 
antisymmetric wavefunctions which are consistent with 
p(r), and the corresponding minimizer characterizes the 
state of "strictly correlated electrons" (SCE). In terms 
of KSDFT, T4e^^[p] is the sum of the Hartree energy 
and the exchange correlation energy, and the exchange- 
correlation functional in the SCE limit can be recovered 

by 



E'.^^[p]-Vl^-[p]- 



1 



p(r)p(r') 



drdr'. (4) 



Eq. (3) allows for treating the kinetic energy func- 
tional and the Coulomb repulsive energy functional on 
the same footing but with different numerical techniques. 
The minimization of the kinetic energy functional Ts[p\ 
gives rise to a energy minimization problem or a nonlin- 
ear eigenvalue problem known as the Kohn-Sham equa- 
tions. Their efficient treatment has been extensively ex- 
plored in the past few decades. The minimization of the 
SCE functional V^^^[p] gives rise to an optimal trans- 
port problem, for which numerical methods are still very 



sparse. From the mathematical point of view, the re- 
sults for the optimal transport problem are primarily con- 
cerned with quadratic cost functions. Proper treatment 
of the Coulomb cost function only appeared recently^^ 
for N = 2. However, a rigorous mathematical general- 
ization for systems with 3 or more electrons has not been 
achieved yet. 

Formally, the optimal transport problem is solved by 
minimizing over all antisymmetric 3A^-dimensional wave- 
functions that are consistent with the given electron den- 
sity p(r). Following physical intuition^^, the optimal 
transport problem can be solved by finding N co-motion 
functions {fi(r), f2(r), . . . , fAr{r)}, f, : ^ E^. Each 
fi(r) represents the optimal position of the z-th elec- 
tron given the position of the first electron at position 
r, with the natural definition that fi(r) = r. Since 
the electrons are indistinguishable and distributed ac- 
cording to the same density p(r), the co- motion func- 
tions should satisfy the mass conservation constraint that 
p{ii{r)) dfi(r) — p{r) dr. Then Vgg™[p] is given in terms 
of the co-motion functions by-'^^'^"^ 



KT^ip] 



N N 



1 



^^|f.(r)-f,(r)| 



dr. 



(5) 



The co-motion functions are implicit functionals of the 
electron density, and can be obtained via semi-analytic 
formulations for spherical symmetric atoms^^'^^ and 
strictly ID systems^"'. However, these semi-analytic for- 
mulations are system specific, and the co-motion func- 
tions cannot be obtained in practice even for general sys- 
tems as simple as a dimer in 3D. 

As an alternative to the co-motion framework, the 
Kantorovich dual formulation of the optimal transport 
problem^'^'^^'^^"^'^ introduces an auxiliary quantity called 
the Kantorovich potential u(r), in which Vg|'~^^[p] can be 
obtained according to 

Kf'W = ^max J u(s)p(s)ds, 

N N N 

s.t. 5:u(r,)<EEur— H' ^^'■'^-1- 

i=l i=l j>i ' * ' 

The Kantorovich dual problem (6) is a linear program- 
ming problem with respect to u, and has the potential of 
treating general systems with an arbitrary electron den- 
sity. However, the Kantorovich problem introduces an 
infinite number of linear constraints due to the arbitrary 
choice of {ri}fLi, and cannot be directly implemented in 
practice. 

Our novel method to overcome the difficulty of infinite 
number of constraints in the Kantorovich problem reads 
as follows. First note that the long-range asymptotic 
behavior of the Kantorovich potential is 



w(r) = ^(r) + C, 



(7) 
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where C is a constant chosen such that the function ^(r) 
vanishes at infinity and satisfies^'^ 



v{r) 



TV- 1 



as r — > oo. 



(8) 



Without loss of generahty we refer to v{r) also as the 
Kantorovich potential in the following discussion. We 
introduce a functional g[v] of v{r) by 



9['"\ 



N N ^ N 

min > > ■; r — > virA 

1=1 j>i I * J I i=i 



(9) 



where the minimization is performed over all possible 
choices of the positions of the N electrons. The Kan- 
torovich dual problem (6) can then be written as 



Vef'^ip] = ^ majc J v{s)p{s) ds + C, 
s.t. g[v] > NC. 



(10) 



Here we have used the normalization condition of the 
electron density, / p{v) dr = N. Eq. (10) is a con- 
strained optimization problem with one inequality con- 
straint. The Karush-Kuhn- Tucker (KKT) condition^ ^ 
implies that the optimal v{y) should satisfy the equal- 
ity constraint 



g[v] = NC. 



(11) 



From Eq. (11) it is straightforward that the constrained 
optimization problem (10) can be converted to a nested 
unconstrained optimization problem by eliminating the 
parameter C, resulting in 



SCEr 



— max 
N V 



v{s)p{s) ds + g[v 



(12) 



We also remark that Eq. (12) can be viewed as a saddle 
point problem 



T/^,^^[p] =maxmin/[i;,{rj] 

^ {ri} 



(13) 



with 



nv,{r^}]^ 



-. » N N ^ N 

-(/.(s)p(s)ds + ^^^^^— ^-^.(r,)). (14) 



Since v{r) has an infinite number of degrees of freedom 
and {r^} has 3A^ degrees of freedom, the saddle point 
problem (13) finds a rank-3iV saddle point of the func- 
tional /[w,{ri}]. To the extent of our knowledge, there 
is no efficient numerical method for finding such a sad- 
dle point directly in such a high dimensional space. In 
practice we solve Eq. (12) via a nested unconstrained op- 
timization approach. 



The minimization (9) to calculate g[v] poses some hid- 
den difficulties: at any set of minimizers {r^}, it must 
hold that 



N 



|ri 



Vt;(ri) = 0. 



(15) 



This is precisely Eq. (20) in Ref. 12. Thus, at the exact 
dual potential w(r), one recovers the co-motion functions 
fi(ri) by fixing ri and minimizing (9) with respect to 
r2,...,rAr. In particular, since (15) holds for arbitrary 
ri, the minimizing set {r^} is not unique. As a conse- 
quence, the functional derivative ^||^(r) cannot be an- 
alytically computed for the exact Kantorovich dual po- 
tential v{r). Thus we use derivative-free methods^^'^^ 
to solve the outer optimization of the Kantorovich dual 
problem (12). The inner optimization (9) and (15) for 
calculating g[v] is a standard optimization problem and 
is solved by the quasi-Newton method. Our numerical re- 
sults indicate that this hybrid approach can indeed solve 
atoms and small molecules with reasonable parameteri- 
zation of the Kantorovich potential. 

Special care should be taken when parameterizing v{r) 
numerically. The long range asymptotic behavior (8) in- 
dicates that the size of the computational domain needed 
to represent v{r) is much larger than the size of the do- 
main to represent the electron density p(r). For smooth 
w(r), it turns out that the correct asymptotic behavior of 
v{r) can be efficiently preserved by introducing a pseu- 
docharge associated to v(r), denoted by m(r), i.e., 



w(r) 



m(r') 



dr'. 



(16) 



The asymptotic behavior (8) translates to the following 
constraint on r7i(r): 



m(r) dr = — 1. 



(17) 



Compared to ^(r) which decays as {N — l)/|r| for large 
|r|, physical intuition suggests that the support size of 
m(r) should be much smaller and is comparable to the 
support size for the electron density p(r), as shall be con- 
firmed by our numerical results below. We remark that 
the co-motion functions are discontinuous for strictly ID 



systems 



13,15 



In particular, Eq. (15) implies that Vw(r) 



is discontinuous for these systems, and the pseudocharge 
will consist of (5-functions and is difficult to discretize. 
Therefore we discretize t;(r) directly on a grid which 
matches the asymptotic condition (8) for strictly ID sys- 
tems. 



III. NUMERICAL RESULTS 

a. Beryllium atom. To illustrate the performance 
of the nested optimization method in practice, we first 



4 



study the beryllium atom with 4 electrons. Similar 
to Ref. 12, the electron density is provided non-self- 
consistently by a configuration interaction calculation 
with Slater-type orbitals^"*'^^. Specifically, p(r) is a 
linear combination of terms r^e~'^^ with j — 0,1,2. 
Since p(r) for beryllium is spherically symmetric, the 
co-motion functions can be obtained semi-analytically^^ 
with numerical optimization performed on the angular 
part of each co-motion function. Our calculation gives 
V;|CE[^] ^ 4 X 0.812132. 

For the Kantorovich dual formulation, we initially 
parametrize the pseudocharge m{r) by a single Gaussian 
function as 



m(r; a) 



TV- 1 



x3/2 



(18) 



with the value a left to be determined in the optimization 
procedure. The corresponding Kantorovich potential has 
the analytic form 



u(r; cr) 



m(r') _ ^ 



1 



erf 



/2 a 



(19) 



The nested optimization method gives V'gg*-^^'^[p] = 4 x 
0.647 with a = 0.8630, and the relative error of V^^^^ is 
20.3%. The result can be significantly improved by pa- 
rameterizing the pseudocharge m by a sum of two con- 
centric Gaussian functions: 



m(r) = (iV-1) cos^{d)- 



(20) 

which yields Vj,^^'^[p] = 4 x 0.7995 with parameters 
CTi = 0.4507, 0-2 = 1-862, d = 0.6872. The relative 
error of Vj.'-'^[p] is significantly reduced to 1.6%, which 
is quite small given that only 3 parameters are employed. 
The corresponding Kantorovich potential v{r) is shown 
in Fig. 1, in comparison to the (numerically) exact po- 
tential obtained via the co-motion formulation. As for 
Ke^^: the potential v{r) with pseudocharge (20) agrees 
remarkably well with the exact potential. 

b. Quantum wire. Next we study a model quantum 
wire system in ID, for which the co-motion formulation 
can also be solved semi- analytically^^. The system con- 
sists of iV = 4 electrons and the Hamiltonian reads 



I N g2 N N N 

(21) 

where Vcxt{x) — ^yuP'x^ is a confining potential and 



2b "^1462 



exp I I cric 



26 



(22) 



is the effective Coulomb interaction. By increasing the 
length scale L = 2lo~^/'^^ the system approaches the SCE 
limit due to the long-range effective Coulomb interaction 




FIG. 1. (Color online) Kantorovich potential v{r) for the 
beryllium atom: co-motion formulation (thick blue solid line) , 
Kantorovich dual formulation with the pseudocharge m(r) 
parametrized by a single Gaussian (thin magenta solid line) 
and by the sum of two Gaussians (red dashed line) . The green 
dot-dashed line shows the asymptotic expansion (8). 



Wb{x). Concretely, as L increases from 4.5 to 14, the 
quantum wire system transforms from a weakly corre- 
lated system with 2 peaks in the electron density to a 
strongly correlated system with 4 peaks in the electron 
density^'', which cannot be described by the local density 
approximation (LDA)'^ of the KS exchange-correlation 
functional. 

We discretize the Hamiltonian by Hermite functions. 
The electron density is represented numerically on a grid 
and is obtained via self-consistent field iterations (SCF). 
In the Kantorovich dual formulation, we avoid the pseu- 
docharge formulation for this example since the deriva- 
tive v' {x) of the exact dual potential is not continu- 
Q^gi3,i5 ^^1^ ^Yi^ pseudocharge consists of ^-functions. In- 
stead, we discretize v{x) directly on a uniform grid. The 
number of grid points is a compromise between accurate 
parametrization of v{x) and feasibility of the optimiza- 
tion (12). We focus on the cases L = 6 and L = 14, 
and choose the grid spacing Axl somewhat heuristically 
as Axg = I and Axu — 4. We allow v{x) at the grid 
points -Ml, -Ml + Axl, - ■ ■ ,Ml with Mg = 7.5 and 
Mi4 — 28 to be determined by the optimization proce- 
dure, and fix v{x) by the asymptotic formula (8) at grid 
points \x\ > Ml- Between grid points, we use piecewise 
cubic Hermite interpolation. Additionally, due to the 
even symmetry v{x) = v{—x) it suffices to optimize v{x) 
for a; > only. The interval [—Ml, Ml] (almost) covers 
the support of the electron density p{x) and corresponds 
to the characteristic shape of p{x), which we try to re- 
produce by the SCF iteration. As starting point for the 
SCF iteration, we convolve the exact v{x) and p{x) from 
the co-motion formulation with a Gaussian with variance 



and J , respectively. We use linear mixing with param- 
eter A = 0.1. 

Fig. 2 shows the Kantorovich potential v{x) and the 
density p(x) obtained via our method (after 15 and 25 
SCF iterations for L = 6 and L — 14, respectively), 
in comparison to the (numerically exact) co-motion for- 
mulation. The Kantorovich dual formulation correctly 
reproduces the strong interaction limit (L = 14) with 4 
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peaks in the electron density. While the results match 



v(x) X L/2 p(x) X L/2 




2x/L 



-6 -4 -2 2 4 6 -4 -2 2 4 



(a) (b) 

FIG. 2. (Color online) Comparison of the Kantorovich po- 
tential v{x) (a) and density p{x) (b) obtained via the "exact" 
co-motion formulation^* (solid lines) and our dual formula- 
tion (dashed lines), respectively. The green (upper) curves 
correspond to L = 6, and the black (lower) curves correspond 
to L = 14. 

quite well for L = 14, one notices a deviation of the den- 
sity p{x) from the co-motion reference for L = 6. This 
observation is also reflected by the values of V^J^^ (after 
the SCF iteration) shown in Table I. Namely, the rela- 



L 


6 14 


"exact" V^^^^ 


1.025 0.3408 


dual-K Ke'^^ 
relative error 


0.9394 0.3381 
8.4% 0.8% 



TABLE I. Vil of the model quantum wire system in ID, for 
the co-motion formulation (reference) and the Kantorovich 
dual formulation. 

tive error of V^^^ for L = 14 is much smaller than for 
L — 6. The deviation is likely due to numerical difficul- 
ties in the maximization (12). As mentioned above, we 
make use of the Nelder-Mead simplex algorithm^'^ which 
is a derivative- free optimization method for the outer op- 
timization. Unfortunately, the results shown in Fig. 2 
depend quite sensitively on the parametrization of v{x), 
e.g., the choices of the above Axj^ and Af^. For differ- 
ent choices, v{x) might acquire local maxima during the 
SCF iteration. Thus further improvements of the opti- 
mization (12) are required, which we leave as work for 
the future. 

c. Trimer molecule. Finally we apply our method 
to a model trimer in 3D, for which the optimal trans- 
port problem cannot be solved with known techniques 
using the co-motion formulation. For simplicity, the elec- 
tron density p{r) is given non-self-consistently by a sum 
of three Gaussian profiles. The normalization of p(r) is 
fixed by the number of electrons N = 2,3,4,5,6. An 
isosurface of the electron density is shown in Fig. 3a. 
We parametrize the pseudocharge by a sum of 3 Gaus- 
sian functions with the same variance a. The centers 
of the Gaussian functions are located on the black lines 
in Fig. 3a connecting — > 1, — > 2 and —> 3 re- 
spectively, with equal distance R from the origin. The 
variance a and the distance R are to be determined by 



m 




(a) (b) 



FIG. 3. (Color online) (a) An isosurface of the electron den- 
sity p(r) (with normalization 1) of a model trimer. p(r) is 
the sum of 3 Gaussian functions centered at the points 1, 2, 3, 
respectively. Each Gaussian has variance |, and each point 
1,2,3 has distance 1 from the origin, (b) Optimized pseu- 
docharge m of the trimer molecule (solid blue) and density 
p(r) as is in (a) (dashed red), plotted along the line connect- 
ing 1 and in (a). The values 0, 1 on the x-axis match the 
corresponding points in (a). 



the optimization procedure. The results are summarized 
in Table II, including Vj,^^. Fig. 3b shows the optimized 



N 


2 


3 


4 


5 


6 


T/SCB 

a 
R 


0.1973 
1.1073 
0.2260 


0.4617 
0.9804 
0.5322 


0.7584 
0.8313 
0.8538 


1.0711 
0.7741 
0.9062 


1.3959 
0.7315 
0.9417 



TABLE II. of the trimer molecule, and corresponding 

optimized pseudocharge parameters a and R. 

pseudocharge m{r) plotted along the line — )■ 1 as in 
Fig. 3a. Along with increasing N, the magnitude of the 
pseudocharge increases as required by the normalization 
condition (17). The shape of the pseudocharge devel- 
ops from a unimodal function for = 2 to a bimodal 
function for = 6, indicating growing influence of the 
distance R. For large N the bimodal pseudocharge is 
biased towards the negative axis around atoms 2 and 3 
where the electron density is larger. The increase of R is 
accompanied by a decrease of cr, and the support size of 
the pseudocharge remains approximately the same as N 
increases, and is comparable to the support size of p{r). 



IV. CONCLUSION AND FUTURE WORK 

In this paper we present a nested optimization method 
for solving the Kantorovich dual problem to obtain the 
exchange correlation functional in the SCE limit for 
strongly correlated systems. With reasonable parame- 
terization which preserves the asymptotic property of the 
Kantorovich potential, the Kantorovich dual solution can 
be obtained for atoms and small molecules. Based on the 
Kantorovich dual formulation, one can combine the SCE 



6 



exchange-correlation functional with existing exchange- 
correlation functionals for the uniform electron gas in or- 
der to improve the performance of KSDFT for strongly 
correlated systems. 

Due to the difficulty in obtaining the functional deriva- 
tive ^^j^ (r) by an analytic formula, in practice the outer 
optimization of the nested optimization method is solved 
by derivative-free optimization methods. However, our 
numerical results indicate that the derivative-free meth- 
ods may get stuck at local minima. Moreover, the 
derivative-free methods are not suitable for optimizing 
with respect to a large number of degrees of freedom. 
More efficient numerical methods need to be developed 



in order to obtain the Kantorovich dual solution for more 
general systems in the future. 
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